{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# One-way ANOVA\n",
    "\n",
    "For the one-way ANOVA, you can use either the function in \"scipy.stats\". Alternatively, you can use the \"statsmodel\" tools: they provide more, and nicer formatted, information.\n",
    "\n",
    "Author:  Thomas Haslwanter, Date:    April-2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "import urllib\n",
    "from statsmodels.formula.api import ols\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The importance of the variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy/EUOrgAAAgAElEQVR4nO2da5Ac13Xf/2dmgeVDpGQvH8WQ2IByyUmpAoOPFZUNTXsZsGASUqSo+EEuR1mYtLEGQrDEUqkkbcmwN0aRiOJEXocADSwkwtjICj8IoqSwKJqP4piowtAMIAAEKUo0LcELGrIhrku0LJoL7szJh57G9jb6cbv79kx3z/9XtTWP7uk+d3r2f0+fe+65oqoghBBSfmq9NoAQQogdKOiEEFIRKOiEEFIRKOiEEFIRKOiEEFIRBnp14ssuu0xXr17dq9MTQkgpOXLkyBuqennQtp4J+urVq3H48OFenZ4QQkqJiPxN2DaGXAghpCJQ0AkhpCJQ0AkhpCJQ0AkhpCJQ0AkhpCJQ0AkhpCJQ0AmpIs0msGOH80j6hp7loRNCcqLZBNatA86eBVauBJ55Bhgd7bVVpAsYeegiclJETojIMREJnA0kImOd7S+LyF/YNZOQPiStl91oOGLeajmPjUYe1pECksRDv1VV3wjaICLvAfAQgNtVdU5ErrBiHSH9ShYve2zM+Yz72bGxPC0lBcJWyOU3AHxdVecAQFXPWDouIf1JkJdtKuijo04H0Gg4Ys5wS99gKugK4EkRUQB7VHXGt/0XAawQkQaASwD8sarO+g8iIhMAJgBgeHg4tdGEVJ6sXvboKIW8DzEV9JtV9XQnlPKUiHxPVZ/zHedGAOsAXAigKSLPq+qr3oN0OoIZABgZGeFipoSEQS+bpMBI0FX1dOfxjIg8CuAmAF5Bfx3AG6r6MwA/E5HnAKwF8Op5ByOEmEEvmyQkNstFRC4WkUvc5wDWA3jJt9s3AdwiIgMichGADwJ4xbaxhBBCwjHx0K8E8KiIuPt/VVWfEJHNAKCqu1X1FRF5AsCLANoAvqSqftEnhBCSI6Lam1D2yMiIcoELQghJhogcUdWRoG2c+k8IIRWBgk4IIRWBgk4IIRWBgk4IIRWBgk4IIRWBgk5I0WAtc5IS1kMnpEiwljnJAD10QooEa5mTDFDQCSkSbpXFep21zEliGHIhpEiwyiLJAAWdkKLBKoskJQy5EEJIRaCgE0JIRaCgE1J0mJdODGEMndin2eSgni2Yl04SQEEndqEA2SUoL53fJwmBIRdilyJPjClj6IJ56SQB9NCJXVwBcj30oghQWe8cmJdOEkBBJ3YpqgCVOXTBvHRiCAWd2KeXAhQ2IBt051Dkwdsi20YKCwWdVIeZGWDrVscLHxxcHlZx7xxmZ53XJ04A992XfwgmjTAnCQ9R+IkHCjqpBs0mcM89wOKi83phITissn+/I5QiQLvt/OUVgkkbtzcND5V1XIDkBrNcSDVoNBxxdqnXzx+Q9Qplu+3skyZ7xDRbJm3Gj2lmS5EzikhPMPLQReQkgJ8CaAFYVNWRkP0+AOB5AB9X1a/ZMpL0KUnCCWNjTphlYQGo1YCdO8//jD+OPj0NzM/nFw5Jm/FjOrBc1Iwi0jOShFxuVdU3wjaKSB3AFwD8eWarCEkaTjARQRsZOF6v+O23nZh82HGynM9kYLmoGUWkZ9iMod8L4ACAD1g8JulXkqYZmnrzWTNwxsacUEirBagC+/YB4+PRop6n0DKlkXgwjaErgCdF5IiITPg3isjVAD4GYLdN40gfYxpHbjaBLVuAW28Ftm1zvPo8Z4KOjgJ33+0MqgLOICxj16QgmHroN6vqaRG5AsBTIvI9VX3Os30awGdVtSXuDz2ATmcwAQDDw8NpbSZFxtRTjtvPJJzghmXeftvxloHuTBoaH1/Klili7JqpjH2LkaCr6unO4xkReRTATQC8gj4C4JGOmF8GYIOILKrqN3zHmQEwAwAjIyOa3XxSKEzj3qb7xYUT3LCMK+Yi3RHYIseumcrY18SGXETkYhG5xH0OYD2Al7z7qOq1qrpaVVcD+BqA/+IXc9IHmKbR2Uq384dlfud37ApYVHri6CgwOdl9sYxLmWQqY19j4qFfCeDRjvc9AOCrqvqEiGwGAFVl3Jw4mKbR2Uq3y9NTLqKna2ITUxn7mlhBV9UfAFgb8H6gkKvqb2Y3i5QSU4E1jY93I2slDNvFvOLaY9JeE5uKHA4iucOp/8QupgIbtZ+JJ5r3wJ9NTzeoPcCS/YDZ3UCUTf7vg0Lel1DQSb6kEd44T7Qb4RCbnq6/PbOzy7NkNm5c2r6wAExNOX+m3ncRw0OkJ1DQSX5EVT+MIs477lZtc1uerr89wHL7Aef9hQWnxszTTwMHDwZ/X0E2lbnWO7EKBZ3kg0n1wzDvPc47LurAn7c9wHL7ve0Blnvo4+PO39SUI+ZJK0AW9fsgXYeCTvIhrvphXJggyjsu4sCftz31upMTv7i4vG1eO4Psn5pyPPO8inmRykNBJ/kQV/0wbZjA6wVPTuZieir8pXkBZ8JTVDZK0gyVqPEIDoQSUNBJXiQNmwwNORNmojzMpIN/YQKYR4aMtz1+Dz1JCCRMmDnwSQygoJP8iBKnRmOpHvnQkNlycEm8+jABtC2M3s7BHye32Wlw4JMYQEEn3SVIUE3FamjICd+oxnu+Ycf0vh9XzzxNW7xhIJuCm6TtpG/hEnSkO7g1SGZnzxdaV6xqtXCxajYdL77Vcvabno4WzLDyu249c2CpnnlQXRSTZea6VTclqO2A2TJ4pK+gh07sEBWX9nqybnllV7x/8hPg93/fiTfX6+FC7Ypnu+0cY34+2p6wGL5bz3zPHkfQ3XrmaSYudStd0N/2o0fNQlSk76Cgk+zECaDXk3UZGADuvRf44heXctXb7XChTiOeYTH8uHrmpiGgoE4j7wFX/8SkrGEjUiko6CQ7cQLoCpJ3IQpV4Nix6Fx1Pxs3Oo9RS76ZYHPikrfTSDPgatIBBE1Mevhh82XwSN9AQSfZiRNAV5BmZx3xcdP57rzTmUgTlqvu4hfK8fHsNieZuATEp1QC6dZBNe0A/PbGhY1IX0JBJ9kxmanoCtL4+PL91qyxUzbWNq69UaLr966ThoWCinaZhmuKvgwe6QkUdGKHtGVzTT4XJ5R5ltIN60zChN69EzGxzz8ZyXv3Eheu4XR/EgAFnRSfKPHKewZlWGcSddfges779ztZO2EZKd52zc0Be/cmuwvhdH/ig4JOykGYeOUdjgnrTEyF/sCBaPu8oR2GUEhGKOik3OSdC560xK/fHnfgN86+LCGUvFdvIqVB1E0j6zIjIyN6+PDhnpybVAy/oKUVuKDjpAnn2LLH1GYW7eorROSIqo4EbaOHnhA6QwUkKBfcTYXctQuYmIg/RtDn5ufThXPSDPxG2RX1g2PRLuKhrwQ9qxjTGSoBjcbSUm7ttrME3po18Rcq6HM7dy7PQpmbc5bVm5/vTo9u8oPjakXEQ98Iug0xpjNUAsbGHA/bnYHaapldqKDPzc8vnxA1M+Nsr9WSrZGaFpMfXNpJUKSS9E21RRuF8cIK+BFDTCoYZj3e6KgTLlmxwhHegQHghReALVuiz+v/3ODgkigODzv54a7Ye9f8zKNdLqY/uNHRpbK969YB27Y5j6zE2H+oauwfgJMATgA4BuBwwPb/BODFzt8hAGvjjnnjjTdqNzl0SPXCC1Xrdefx0KH0x3nggeSfD/pc2mOVAn/jbF0A7/GjjnfokOrmzaorVqg6E+RVBwfjzxt2oS68UFXEOU6ttnTOKDtML3DUfkl+JA884NgBOI8PPBD/GVI6gjTY/Usi6JdFbP93AH6u8/wOAH8Zd8xuC7pq7wQ06H/etr4ViqDG2RYbk+M98MCSCAPO87Tn3bPH6RxEVAcGnNdRdphcYLfTGRy080Oo9I+KuEQJupUYuqoe8rx8HsA1No5rm15NrAsL91Q2Hh/UYNuDdybHGxtzQihnzzqvs5x3ft4Jtbjdg1vmN81MUmBpUMdbgTLrD4HlAPoeU0FXAE+KiALYo6ozEfv+FoBvB20QkQkAEwAwPDycxM5SE/Y/X9nkhKAGm4qNaQ63aUGwRmOptkqWErNhFzHpTFIXV/BdMRex80NgOYC+xmhikYj8C1U9LSJXAHgKwL2q+lzAfrcCeAjAL6tq5JIy/TaxKEiXKp3TnqZx/lSkqDoovSBpm0xXcarXnXK4rGlODIiaWJR4pqiITAH4J1X9H773fwnAowDuUNVX447Tb4Le15gK4Y4dToZGq+WI3Lp1joi7r7dvX74Ic9mpdI9O8iLTTFERuRhATVV/2nm+HsAf+PYZBvB1AP/ZRMyzUuT/gyLb1hOSTABwwxQLC04I4rrrzOqg2LCxFxfNdniEP76+xySGfiWAR8VZ3HcAwFdV9QkR2QwAqrobwO8BGALwUGe/xbAeJCtFnq3ZTdsK+b8bZFSS2Vijo06YZetWZ/8HH3ReHz2ar81F/UGFERa/K1s7iHViBV1VfwBgbcD7uz3PfxvAb9s1LZgiz9bslm2F/N8NMyppdoubTeJO3jl6dHl98bSNDesBi/yDCiLsey5bO0gulG6maJFna3bLNhuzXq0TZpSbBbJ9u5kY+79EIHtjXREMmkHZ7R9Us+nMWo2buRpG2Pdc5H8M0jVKV8ulyKm23bKtkPWYooxKEisOqk2SdeGHKO/V1kUziYE1m852Ny9+3z7g2WeTnTNp+iTpK1gPvaSUJoZehOPmHaMyPf6OHcDnP7889/z++5Nn7hTy4pNuwXroKSny/00h54/kZVTW4+btvZrGr23NXC3kxSdFgIIeQiEHHi1T5A4rkjSG5ymCpjEwmzNXCQmAgu7D1Yq5ufRJA0n0ptui6p5vaKhYkzCNsd3T2rgASe4A6F2THOkbQTcds/LOxh7ofDtJ7oyj9MbWkpVp8Z5PZHl2YGmy3Gym59m8ABRqUgD6QtBN/2+9WgEAmzY5axskcd7C9CbIBu++CwvA1JTzZ1sXgu46ajWn07JVE6pr2Ezx6dYFIKRL9IWgJxmz8mpFXIgzyOtPUk3VO9O93QaeftqZ6W7TU4+665ieTrY8ZiFi7jYHOMMuQNIvhpCiEFYoPe+/bi5wkaTuv+niMf5j7tmzfFvYojd+Gw4dUl2/3lkEx10MZ/16e2sT+Ndf2Lw5/YpLlVw7IegCrFhRjYZWekms/gVZVyzK4y+roCf9rWb9bfsFbfPmJaGM0wD33F7RDzq2u7iOSHotyWvlt0qvbub9kgYGlsS9zA2tbA9MogS9lCGXNGNZWces/CETYCm0IuK8HzTAaGKrW5Nqy5alBXEWFpKP94WdKy5CYbKmRCFnp9rC+yX503/K2lDWdulLSinoUb/VvOK8QfH18fF4DTD9v5r3LQdSryfXkrBzRXVm/hj7hg3At7/tLHLv7RQqP7Pc+yWtWVP+hla6ByZhlFLQw36reacBbtzoPF5//dL/uztrO0wDTP+vxsaAwUHHM6/VgJ07k9sedK64Ds7bCbRawDe+sbQtqOyJ9xiFGCTNgyqkIFa+ByZBlFLQw36red1lzswslegeGHBCLF4P1j130P9NmK1+MbTx/xdU1yqug3M7Ae9axUB8OmM/zKQtPVXomEgywoLref/lkeWSdRwoaOB0z56lMTJ3wNIdvHSzRpKe0+Z4VdRgr+lA5qFDTjsGB539Vq50Xkdl+ngHhcs8dkhI2UDVBkXDyOLlBnmcAHDPPc5gp0ut5njp77zjeLF/93fJ7wps3UnEeclJSoyMji6NCUQNnmadSUsIyY9KCTqQPs4bJLLA+WL+0EPOczcE8/jjyYXN1nhVXMeQtIOLu0O3MZOWEJIflRN0L2nWJ/aLrHegctcuYGLCKWvt1kFptZILm63xKpOOIU6kkwxsJp1JSwjpLpUW9CShjTCRDXrPhrDZGK/K2jEkHdhk4gQhxaaygt5sOsWokoRDgkQ27L2iCFuWjiFJPr/3ddIFdggh3aGSgu4fvNu0yX54oAoZYab5/NPTJa2dTkifUUlB9w/eDQ/nP5O0jJjm8x84wFnkhJSBmslOInJSRE6IyDEROW9lZ3H4XyLymoi8KCI32DfVnKEhZxCzVgv2PLdtcx6bzXztaDadAdS8z5OF0dGlEIprq+u51+vO4513Ln89NFT8dhHSjyTx0G9V1TdCtt0B4H2dvw8C+JPOY64EedvNphMecBdxmJ7OfyZpmG1ZZ1J2624iyFa/5+6WNijt0nWE9AFGHroBHwUw25nI9DyA94jIVZaOHUiYt+2KdrvtzO30Fr3ye555epphee2mzMwAv/qrwO/+rtndRJa7gbCObnJyeR2XyUnn+8zSLkJIfph66ArgSRFRAHtUdca3/WoApzyvX++896PsJgYTJkJRudlRVVJte5pZJg/NzDildN1JTXGldLPeDcTZ6r1TYBE/QoqLqaDfrKqnReQKAE+JyPdU9TnPdgn4jPrfEJEJABMAMDw8nNhYL2HCEpdS6Gan7NiRb/glbWpjs3l+uYG4UrpZQ0lRtpqEYwghxcBI0FX1dOfxjIg8CuAmAF5Bfx3AKs/rawCcDjjODIAZABgZGTlP8JMQJUImKYXd8DTTpDY2GueXG3BL6YbF1E3bEhWTD7M1qLPwhmLijksI6R6xgi4iFwOoqepPO8/XA/gD327fArBVRB6BMxj6pqrmFm5x8Q52el+bfraInqa/LrpbbiAqrGLSlrRhGZNwDMvoElIMTDz0KwE8KiLu/l9V1SdEZDMAqOpuAI8D2ADgNQBvAbgrH3OXk1VMXK/UHVBMO33eZqcQJM7NJjA1tbQ4/dtvA7Ozzv7e/ZIUH5udDQ+x+I8Z1Vl0M3OIEBJDWF3dvP9s1EM3XdE+qmZ4ltrk3ViH1z2HtyY74CxK7dYv95/brW/urWnutXXlyuDPpmkP1yImpLugqvXQveGAeh3Yt+/8lYRmZ89/3+tBZvEwu+GdetMwRZZWFVpcdB5Vl5/bnRjkLmS9bx/w7LPLPe25OWDv3vPtTtOeooauCOlHSi3oUSI1Owvs3798abUgkcoyONqNgVV31qvq8uXv6vXlS+G55240nMU3XIIWi242ne/Gb3fa9lShrg0hVaDUgg6EixTgPHfFPGyNzCweZt7eqX/W64MPLl+MGggu7btixZKHnqTN9LYJKTeimil7MDUjIyN6+PB5ZWEy4R3QA5ZXXLz7bvOKi2kHOm0PkO7Y4cyEbbWcNmzfHl261j3/0BBw9KjzHhehIKRaiMgRVR0J2lZ6D92L/9Y/7cSeNJkzeaTvJQmBMH2QEFIpQfeTdmJPmoHOPAZIk4RAmD5ICKm0oJvgD5OkHRjMa4DUpFNKszpTVjg7lJDiUTpBtykkYWGKNKGaXg0odmN1pqhzMrxDSHEolaDbFpKwMEXaNLxepO952wAsX53Ji82OkOEdEkTzVBONkw2MrR7D6Cr+IHpBqQTdtpBUoRSsSRuydIRBHYF/QtfcnLMfRb1/aZ5qYt3sOpxtncXK+ko8M/4MRb0HlErQbQtwUfKus3jPcaVv3UlXaTrCsI7APac7C3fvXmcOAEMv/UvjZANnW2fR0hbOts6icbJBQe8BpRL0PAS417McbYSRgtrgj62nGTCNuiNySwUsLjL0QoCx1WNYWV95zkMfWz3Wa5P6klIJOmBPgIuSpZFXPNofW9+0yYmvJ2lv3B1RFUJW/UZece7RVaN4ZvwZxtB7TOkE3QYzM8DWrY7YDQ72NlSQRhRNOiP/cdNkvsTdERUlZEXMyDvOPbpqtO+FvNcDw30n6O4Sb261wrj1Om2cL0rwkoqiaYjGhtiadBy9DlkRcxjnzpciDAz3naA3GsHrdeYRgkkivnmU7M0itsw1rx5Z4ty99jzzwHabitBh9p2g+5d427nTeT8P8cojPt6tuDVzzatH2jh3ETxP2+TRpiIMDPedoAeFInbsyEe88hDfbsWtOeBZTdLEufP0PHvl+efRpiIMDPedoAPnhyLyrMOSh/h2I27NAc9qk0RI8/I8e+n559WmXg8M96Wg+8lTvMo8aFhm20k4SYU0L8+zlzHnInjTeUBB72BbvIqS506InzRCmofn2euYc6+96TyojKAXSUCZIUKKjF9Ihy4awo6DOxIPlGb1bqvqJfeSSgh60QS06BkiRer8SPfxCunQRUO474n7EsWx/SGb6dunMf/WfCpRjvOS03QcVUyxNKUSgp6HgGYRvahB1l6LadE6P9IbXCHdcXBH4vCLN2SzsLiArY9vRVvb1gc20wyaVjHFMgk10x1FpC4iR0XksYBt7xaR/ysix0XkZRG5y66Z0bgCWq/byVJxRW/bNuex2Uz2eXeQdfv25YKZ9bg2COr8SP/ihl/qUke9Vsfcm3Nonor+YXo/U6vV0NLWsg7BFkGxftufaZ5qYsfBHbFtLgtJPPRPAngFwKUB2+4B8F1V/Q8icjmA74vIn6nqWRtGxmE7S8WGxx80yFqEUAzzy4kXN/wye3wW+47tw97v7MX+4/sjPduokI3Ngc00g6ZJPhPmzZc5ZGMk6CJyDYAPAbgfwKcCdlEAl4iIAHgXgH8AsGjLSBNsZqnkJXpFEFPmlxM/o6tG0TjZwGJ7MTb04hW7yVsmAQBrrlhTmAqOST4T5s2XOWRj6qFPA/gMgEtCtu8E8C0Apzv7fFxV2/6dRGQCwAQADA8PJza2W+Q5IagIYsr8cuLHxLMN82jzTP9Lc2zTzwS12WZufC88/VhBF5EPAzijqkdEZCxkt18DcAzAvwfwCwCeEpGDqvqP3p1UdQbADACMjIxoFsPzJi/Rc4/pxq4prKQImHi2RSg+ZZOwNtvIjZ85MoOtj29FS1sYrA92zdM38dBvBvAREdkA4AIAl4rIV1T1E5597gLw31RVAbwmIj8E8K8BvGDd4pLDLBNSVII8W6+X2euJQHngb7ON3PjmqSbuefweLLadqPNCa+FcOCdvjz1W0FV1EsAkAHQ89E/7xBwA5gCsA3BQRK4E8K8A/MCuqdWgCAOjhJgQFGIJE7syDyT6yRpCapxsoO2p0V2XOoYuGupKbD51HrqIbAYAVd0NYDuAPxWREwAEwGdV9Q07JlaLIgyMEmJCUIhl8pbJQC++TAOJeXc+Y6vHMDgwiIXFBdRqNezcsBPzb813JVyVSNBVtQGg0Xm+2/P+aQDrbRpWVYoyMEpIHHEhFlcY596cixWronjw3eh8gsI2zVPNroSrKjFTtGwwy4SUgah4slcY67U6BmoDQBuhwm8qonkLf7cGdvOIzZtAQSekD7BdE8UrjGgDm27YhOF3DwcK/1RjCgutBbS1HZvjnsZ7LkJtdxO6Ud2Rgk5IxcmjJopfGMfXjofG1hcWF9BGGzWpRYpoGu85bW332eOzkcf1nyNr1ku3wk0UdEIqThqhjPtMkrz1NtqooYbbrr0NU2NToedO4z2nDaHsP74fZ1tnY8scRHUYJkLd7QFjCjohFcdWTRS/gMWFEPzHiBJzIF2cOe9OIGxfU6Hu9mQsCjohFcdGTRQgeY2TtOdNInh5dwJh+5oKdbdj9hR0QvoAr1CaxnS9n0lTN91/DFOSxpzz7gQ2rt0IAMvGCUyFOk3MPgsUdEIqjF8c08Z0u+Vp5h1zDqoWaWrL9Vddv+y7TDJr1jRmnxUKOiEVJUgc08Z0u5VHnWfM2aSz8IqxycpMJrNmuxlHp6ATUlGChCSLp92NPOq87gRM8uGD1kp1bQGAxfYiFBopyra/86RQ0AmpKK6QLCwuQEQwdNGQdU/bdo61176hi4bOVSnMcmzTfHi/GM+/NX8u/v3lo192JlEBGKgNhIqyX7zdNmRZSDsJFHRCKsroqlFM3z59ri73fU/chzVXrLHmaecV7x5dNYoTZ05kqiceFDqJy4cP8qRHVzmrObU76/UIBHddd1eoLf4Oybs8XzeKlhkvEk0IKR/zb82jre1lYQYvWRZJDlvCLSl+G9x64u+030Fb28vqiZsca8tjW3Dr/lux7dltWDe7DkMXDZ1b1HpwYDA0H94V4+23bl8mvt5FsS8YuADja8cjbRhdNYrJWyYDKyzmDT10QipMVPw2q4dtIzYcNojoryducmz3WG8vvg2FsyCaN3SSNFXT+16aMFXY95NnKQAKOiEVJkqMomZBmopf1nh82CCiv564ybHdY7liLpBloZMs4ul+3r2bMGlv0PeTd1omBZ2QihMmZmHT+90BxFqthl0bdmHixonExzYlLG6dpqMYumgINalBoRioDeDu6+4OLBqWliAxBqKXlfN/P3mnMFLQCelTgoRzy2Nb8M+L/wwAaLfb2Pr41nMDqd2ywX0/6SDofU/ch1a7hVqthgfveDCyIwr6vN+T9tvkF+PZ47PnJgyZett5pzBS0AnpY/wlAR4+9vCy7S1tYfb4bK4TiuLE2yQE5M1kERXMvzVvfP6g/POg7BS/GANI7G3nPUGLgk4IAeCIYqvdOvdaIBioDWDfsX1YbC+mjvlmGQQ0jTln8Xz9nveB7x4IFOqggmWuh16v1TH35hyap5pGos6ZooSQXPGKYr1Wx93X3Q0A2PudvaljvlkHAeNizt7OIq3n6429r6yvxJ3vvxMH5w4Gdg5+MXYnHu07tg97v7M391otcVDQCSEAwrMyvHHipDHfrIOAUZ73zJGZ8yYfxRXc8uOPvU/fPo2JGyew5oo1xpk+jZMNLLYXu1bzPAoKOiF9gknow++Bpo35uudyJ/Wk7RDCys+6k48W24sAcG7yUVIh9cbe0QYOfPdA4tm0vVyn1A8FnZA+IEvoI03GiX+QMWsdk/3H92NhcQFfOvol7Nqwy5kBm2LykR9vvZs22nj6h0/j4NzBxN9PNypRmmA89V9E6iJyVEQeC9k+JiLHRORlEfkLeyYSQrJia5p+mnPNvzWPyVsmUwtd42TjnOAuthex9fGtGLpoCIMDg6ihhoHaQOjko7jSBq4Y3/be21CTWmiJhDjc6f69FHMgmYf+SQCvALjUv0FE3gPgIQC3q+qciFxhyT5CiAW6GRawfa6x1WOo1WrnPPKWtoym85velYyuGsXU2FToQKhN8pz2DxgKuiD7oM4AAAi6SURBVIhcA+BDAO4H8KmAXX4DwNdVdQ4AVPWMNQsJIZnpZljA9rlGV41i14Zd5wZAB2oDmHtzDgAiB0GTDMiGDQjb/L7ynvYPmHvo0wA+A+CSkO2/CGCFiDQ6+/yxqnZnET1CiBF55j/nfS438yRJimDSOwX/JKuk4hvXAXRj5aJYQReRDwM4o6pHRGQs4jg3AlgH4EIATRF5XlVf9R1rAsAEAAwPD2exmxBSIUwzcMJSBMM+H7TAswlJxdekA+hG2MvEQ78ZwEdEZAOACwBcKiJfUdVPePZ5HcAbqvozAD8TkecArAWwTNBVdQbADACMjIyojQYQQsqHV4ABGBcEiyoo5i+a5b5Xr9XPfT6PVEST5e3cc+cd9ooVdFWdBDAJOJksAD7tE3MA+CaAnSIyAGAlgA8C+CO7phJCqoBfgDeu3XguiyWuIFiQKO44uCMwg8d9r9VqYc+RPYlmcZqKr+nydt7j5hn2Sp2HLiKbAUBVd6vqKyLyBIAXAbQBfElVX7JkIyGkhISFQfzhDADnZbHEDWC6xwHCvemV9ZXnFrtQKBYWFzDVmApdsSjoPHH7mS5v1y0SCbqqNgA0Os93+7b9IYA/tGUYIaS8hMWUm6eamHtzDgO1AaDtiO742nFcf9X1y6bxx4U4/McO8qa9dVbeab2TeuJQFP7OpJdiDnCmKCHEMmExZWB5XHvTDZvODVaOrho1rp8SNGAZNKnHPe742nFMNabw9A+fjoxxp6FIs0QBCjohxCJRMWWvEKMNDL97+Ly6MSaCmCYdMc+JQ91MB42Dgk4IsTaJJi6mbCNtL41XXDRPOi9EtTfZgyMjI3r48OGenJsQsoTNGYxxx8p76nu36GU7ROSIqo4EbaOHTkifY3MGY5wnXKTwRFq6MYU/LRR0Qvoc2zMYqyDaUXRjCn9aKOiE9Dn9El+2RZEWtPDDGDohpCtUJX4OMIZOCOljihx3TkNRw0rGKxYRQqpF3Go+Nunmikn9DD10QvqQbnvMRY47VwkKOiF9SLczNWwNvFYpDp8HFHRC+pBeeMxu3NkN9SQV5arF4fOAgk5IH9KrVMUsolzk/O+iQEEnpE/pRaZGFlFmHD4eCjohpGtkEWVOgIqHE4sIIV2FA5vZ4MQiQkhhKOqknCrAiUWEEFIRKOiEEFIRKOiEkER0s2QASQZj6IQQYzi5p9jQQyeEGMMiW8WGgk4IMcbNI69LnZN7CohxyEVE6gAOA/hbVf1wyD4fAPA8gI+r6tfsmEgIKQqc3FNsksTQPwngFQCXBm3sCP4XAPy5BbsIIQXFVh45JxjZx0jQReQaAB8CcD+AT4Xsdi+AAwA+YMc0QkhV4eBqPpjG0KcBfAZAO2ijiFwN4GMAdkcdREQmROSwiBz+8Y9/nMhQQkh14OBqPsQKuoh8GMAZVT0Ssds0gM+qaivqWKo6o6ojqjpy+eWXJzSVEFIVOLiaDyYhl5sBfERENgC4AMClIvIVVf2EZ58RAI+ICABcBmCDiCyq6jesW0wIKT0cXM2HRNUWRWQMwKfDslw6+/wpgMfislxYbZEQQpITVW0xdR66iGwWkc3pzSKEEGKTRFP/VbUBoNF5HjgAqqq/mdUoQgghyeFMUUIIqQgUdEIIqQgUdEIIqQgUdEIIqQg9WyRaRH4M4G9SfvwyAG9YNKeXsC3FhG0pJmwL8C9VNXBmZs8EPQsicjgsD7NssC3FhG0pJmxLNAy5EEJIRaCgE0JIRSiroM/02gCLsC3FhG0pJmxLBKWMoRNCCDmfsnrohBBCfFDQCSGkIpRO0EXkdhH5voi8JiKf67U9SRGRkyJyQkSOicjhzns/LyJPichfdR5/rtd2BiEiD4vIGRF5yfNeqO0iMtm5Tt8XkV/rjdXBhLRlSkT+tnNtjnXWAHC3FbItIrJKRJ4VkVdE5GUR+WTn/dJdl4i2lPG6XCAiL4jI8U5b/mvn/Xyvi6qW5g9AHcBfA3gvgJUAjgN4f6/tStiGkwAu87333wF8rvP8cwC+0Gs7Q2z/FQA3AHgpznYA7+9cn0EA13auW73XbYhpyxScev/+fQvbFgBXAbih8/wSAK927C3ddYloSxmviwB4V+f5CgB/CeDf5n1dyuah3wTgNVX9gaqeBfAIgI/22CYbfBTA/s7z/QD+Yw9tCUVVnwPwD763w2z/KIBHVHVBVX8I4DU4168QhLQljMK2RVV/pKrf6Tz/KYBXAFyNEl6XiLaEUeS2qKr+U+flis6fIufrUjZBvxrAKc/r1xF9wYuIAnhSRI6IyETnvStV9UeA86MGcEXPrEtOmO1lvVZbReTFTkjGvR0uRVtEZDWA6+F4g6W+Lr62ACW8LiJSF5FjAM4AeEpVc78uZRN0CXivbHmXN6vqDQDuAHCPiPxKrw3KiTJeqz8B8AsArgPwIwD/s/N+4dsiIu8CcADAfar6j1G7BrxX9LaU8rqoaktVrwNwDYCbROTfROxupS1lE/TXAazyvL4GwOke2ZIKVT3deTwD4FE4t1V/LyJXAUDn8UzvLExMmO2lu1aq+vedf8I2gL1YuuUtdFtEZAUcAfwzVf165+1SXpegtpT1urio6k/grPR2O3K+LmUT9P8H4H0icq2IrATw6wC+1WObjBGRi0XkEvc5gPUAXoLTho2d3TYC+GZvLExFmO3fAvDrIjIoItcCeB+AF3pgnzHuP1qHj8G5NkCB2yIiAuDLAF5R1S96NpXuuoS1paTX5XIReU/n+YUAbgPwPeR9XXo9Gpxi9HgDnNHvvwbw+V7bk9D298IZyT4O4GXXfgBDAJ4B8Fedx5/vta0h9v8fOLe878DxKH4rynYAn+9cp+8DuKPX9hu05X8DOAHgxc4/2FVFbwuAX4Zza/4igGOdvw1lvC4RbSnjdfklAEc7Nr8E4Pc67+d6XTj1nxBCKkLZQi6EEEJCoKATQkhFoKATQkhFoKATQkhFoKATQkhFoKATQkhFoKATQkhF+P9GcG8CeCHePgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "centers = [5, 5.3, 4.7]\n",
    "std1 = 0.1\n",
    "colors = 'brg'\n",
    "\n",
    "data1 = []\n",
    "for ii in range(3):\n",
    "    data1.append(stats.norm(centers[ii], std1).rvs(100))\n",
    "    plt.plot(np.arange(len(data1[ii]))+ii*len(data1[0]), data1[ii], '.', color=colors[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy/EUOrgAAAgAElEQVR4nO2df4wex3nfv8/70sdQsZ3EJylxZbKUCyOAIEaqeDB8cK0eSsWlBbVKkbZwgPRY0dVBqJiaRYNWV8PNqVbv+jM4AxYKUo0ZsjHiFGEKG4Gq2jr4KhZ8o+RoS/6lqlZs6iibtZhL66RVdOe79+kf+y65XO7vnd2Z2ff7IQ4v333fd3Zmdvc7zzzzzIyoKgghhPhLz3YGCCGE1INCTgghnkMhJ4QQz6GQE0KI51DICSHEc3bZOOnNN9+s+/fvt3FqQgjxlgsXLvyRqt4SP25FyPfv34+1tTUbpyaEEG8RkVeTjtO1QgghnkMhJ4QQz6GQE0KI51DICSHEcyjkhBDiORRyQgjxHAo5cZfBAFhaCl4JIakUjiMXkU8DeADA66p65+jYOwD8FoD9AC4C+Nuq+r/NZ5OMHYMBcOgQsLUFTEwAKyvA9LTtXBHiJGUs8l8HcDh27DEAK6r6HgAro/eE1Gd1NRDxnZ3gdXXVdo4IcZbCQq6qzwH449jhBwGcHv3/NICfM5QvMu7MzASWeL8fvM7M2M4RIc5Sd4r+T6rqZQBQ1csicmvaF0VkDsAcAOzbt6/maUnnmZ4O3Cmrq4GI061CSCqtrbWiqicBnASAqakp7i9H8pmepoATUoC6USvfF5F3AsDo9fX6WSKEEFKGukL+eQBHRv8/AuBzNdMjhBBSksJCLiK/CWAA4KdF5DUR+QiAfwngZ0XkWwB+dvSeEEJIixT2kavqL6R8dMhQXgghhFSAMzsJIcRzKOSEEOI5FHJCCPEcCjkhZeBCXsRBrGy+TIiXcCEv4ii0yAkpChfyIo5CISd+UcW1YcodwoW8iKPQtUL8oYprw6Q7hAt5EUehkBN/SHJt5Ilpld9kwYW8iIPQtUL8oYprw7Y7hFEupAVokRMMBp54C6q4Nmy6Q7LcOt5UOvEBCvmY411EXRXXhi13SJpbx7tKJ65D18qYw4i6Bklz68Qr/cwZul9ILWiRjzmh1oTGISPqUqjiCklz60Qrvd8HTp0CtrdpnZPKUMjHHEbUFaCOKyTJrROt9PV14KmnzEXVkLGEQk4quZDHaqzOdAhjtPIA4PRpdolILSjkpDRjN1Zn0v+UVHnsEo0Fg0sDrF5cxcz+GUzvNXudKeSkNKYNVOcx6X9Kqrz5+Y5XIBlcGuDQmUPY2tnCRH8CK7MrRsWcQk5KM5YDpKZCGMey8tykSQs5zurFVWztbGFHd7C1s4XVi6sUcmIXDpDWgJXnBE1byHFm9s9goj9x9Xwz+2eMpk8h7whtDz5yyZEasPKs07SFHGd67zRWZlfoI/eJtkV17AYfCalJ0xZyEtN7pxtrLCjkhrEhqmM3+EhITZq2kNuGQm4YG6LK8TNCytOkhdw2RoRcRP4hgL8HQAF8DcBDqvqmibR9w4aocvyMkPGmtpCLyG0A/gGAO1T1z0TkPwH4MIBfr5u2j9gS1brjZ2M1U5O0RpshfuOMKdfKLgB7ROSHAG4C8D1D6XqJb0EJHCwlTdB2iN84U3sZW1X9LoB/C2AdwGUAP1DVL9RNl7QHl7IlTZAU4keaobaQi8hPAHgQwO0A/hyAHxWRX0z43pyIrInI2pUrV+qelhjE9m5opJuEIX596bcW4jeuiKrWS0DkbwE4rKofGb2fBfA+Vf37ab+ZmprStbW1WuclZqGPnDQBfeRmEZELqjoVP27CR74O4H0ichOAPwNwCABVumFMC691v368QEkFtNXasJWrTJUQP4p/eWoLuao+LyK/DeDLALYBfAXAybrpknQ6NzgZL9DyMnD8+PUFBOwUunOV7SaheE/eNInjzxznAGlJjEStqOqvAPgVE2mRfLImHXlpPMYLdPZs8uirjemrnDbbCFGrG8DV6BYRwVCHGOqwlTVQugJndjpEURFOm3TkrfEYL9DddwNf+hKgen0BbUxf5bRZ48TDEo/cdeRqdEtPe+j3+hAIB0hLQCF3hDIinDbpyFvjMVqgycnArbKzA/R6gZslLIStmVacNmuUeFgigOsWsFo+vIyNNzboIy8BhbwiJl0YgwGwsABsbgLDYTERThqc9Np4DAu0tBQUYDgERICNjRu/YytvLuGlDy0gvvLg7F2zmL1rlgOcNaCQV8CkCyNMKxTxXq+6CGcZj40+9yYTt9ka+SCOgwFw5gxw6hSwve2ZDy0gbeVBGwLemQgZVW397+DBg+oy58+rLi4Gr0ksLqr2+6pA8Lq4WP1c0bR6PdUPfjD9vFU5f151z57gPHv2GE6/icTzLkATNFpJhgjzKBLcMCZuwDKnXz+vi88t6vl1B+umAufXz+ueJ/Zo//G+7nlijxflArCmCZpKizxGEWu7qtGYZPDF01pYMGNcRc/VqO+8icRtuDJ8GGAI8xhO4hNprdfSxXVT2t4lqEko5DGKPM9Vxr/SGogmxtKSwrIb81Z47ZiPUKQctl0v0Tz2+8DRo8DsbCt56ZLohdjYJagpKOQxiupSWaMxq4EwbYDGz7Wx0WDgRVeiOvLK4UJsp8W67pLohXRplyAKeQzTz0poxE1Otme4JjVGjXorXIzqqEJWOVxxvViq6y6JXpSu7BJEIU/A1LOS5OLY2GjemOqKkVwa0zGh0bS64kKqQVdEr4tQyEtQVieSXBzz8w1ncoQVw82mD7mJmNB4WmPZOubjWwifb/ktAoW8IFV0oq4RZ3tsrRS2fcgmXR9paXXFhWQQ36JZfMtvUWpvLNEmg0Ew8W8waP/cSc92HqER94lPlNe1UBc//vHg1UaZS1GlgkxicncM7rRRGN92AfItv0XxxiJPMviA9izWstZ11Jqu4k5xZWytMLZ9yCZdH02NeFvsWjXlTvAtmsW3/BYmaZZQ039VZnbGZ1M+8kjzE/HiEwyLTjg0MUnQxYmGueU/fz64MI88kv4lG7M22yCtXA5cyKZnMPo249O3/EaB7zM74wYf0KzFmjWBJw8T1nQdo7AJA7CwC/z06eBLp0/f+CXbfvSmyCqXA12rvMk8da1136JZfMtvEbwR8riwAdc0o4mefJ3nz5SXocrYWlNaWag+8r7kgKg1Qla5bLuckO1O6Org37jhjZADNwpbXTdmluVa5/mzGanWlFYWqo+8LzkgakYoE2Pe1s2QcTNnTebp4tT7cUQ0XICnRaampnRtze7+zEUsVwfGqErTpPeiUH3kfcnHSo2SVsGextDTIk/HxXhzEbmgqlPx415Z5CYpYrn6GDbc5Jrkheoj70s+VmoUF2PMa3TDujr1vi6+NXBjK+RJvWHbRpWpcydpihPjjL5b44Cb7qGaeerS4J8pK9o3l9NYCHmSfiQNnlYVurL6FP9+GyJrfZzRiZbEAC5O1XcxTxYwaUX7Fm/eeSHP0o+o5RpuFVlW6MrqU9L32xBZ64akrZYk2mqG+TAxYajJvFfpufjusjKASSvaN5dT54W8qH5UFbqy+pT0/TZE1rrRZqMlibaa/X6wo06ZfS5tuILq9ly64L6qiGkruozLyfbAaOeFvMxGEVWErqw+pa0V3obIWjXabLQk0VZzOAyOqRZrcW25gur0XLrivqqILSvahYFRI0IuIj8O4D8AuBOAAjiqqs4s83TkSPCatytWFaFL8rUvLaVrVZqeudQzbsyoa7uQ0VYzbpHntbi2XEF1ei7WB0LsEbWI5z8w3+q5s1w6bVnqpizyTwJ4RlX/pohMALjJULq1iBsos7PNnCfUp6IGkUuiHcdbo67oiHZ8lDmtxbI1qFCn52Igz7ZdBFWwbRHHXTqTN01i6dwSJm+axPFnjreSr9pCLiJvB3AvgL8LAKq6BWCrbromaNtA6YJBZLwMbUwQKjqiHb7P+034PVuDClVb+pp5ti2IVbEdKhh16UTFW0Qw1CGGOmw8XyYs8ncDuALglIjcBeACgI+q6v+LfklE5gDMAcC+ffsMnDYfk0ZVEb2xHhligCJlKKy9eWJZxvzPOmmV1qfIb1zuOqVRMM9JlncRQXTRYnchVDAcGF06t3S1DnvaQ7/Xh0Aaz5cJId8F4B4Av6Sqz4vIJwE8BuDj0S+p6kkAJ4Fgir6B8+Ziyqgq4zLxPZw3rwylXC95YllUgPNOWqUFnZwEer1g8NPXVrciaZZ3niC6YrHHGxOXQgXjdbh8eBkbb2x44SN/DcBrqvr86P1vIxByJzBhVJUx+Hw04uJklaGU8ZslsIMBsL4O7BrdgllimnfSsi3oYAAcPx6k1+sFu2L7ftFKkGZ55wliEy6MPAs//nlaY+LK7FRbjUptIVfV/yUil0Tkp1X1ZQCHAHyzftbcoQsuk6LkuU1K1UVUYCcng9eQaHz3ww9nhxQVOWmZFjRsGIbDIJJlY6PY7zpCluWdJYimXRh5Fn7S57b94UWw0aiYilr5JQCfGUWsfBvAQ4bSdYIuuEyKcPIkcOxYYKju3p3sNildF+EXoq6RI0euWdgAsG9ffkJFY0iL4GPLbGBQOGrdVrEaTVubeaKc9HndxsRFH78JjAi5qr4A4IalFbtEF1wmWQwGwKOPBmHWALC5me42KV0XcdcIUFxIm4gh9a1lNhATmmTdVom3TrM2qwhknignfV6nMani4/dF+Ds/s5MUY3X12uRHIPB41DZUQytycvJ64Z6dDf6KCGlTMZ0+tcwG6qBJl0TVQdA8UU77vKrromwduDK4W4ROCXlXl5loo1wzM4E7ZXMzGP/71KcMxI9Hrcjl5cAXHZ+wUyRjvrlBkqhzEQ3UQZMhekkCGR7Ps2TzRNmkv7lsHUTL9eb2mzjz4hkKedO0OSOxyjNZ9Tluq1zGvQ1xK3JjA5jP6coXmZ3pYwtd9yIaqAOT/u24uyFpZqOLlmxeHSSVq9/rY2dnBwrFqRdOYfauWSfKEsdLIU963pvqgSedu+wzWWWp27B8Z84Ab75ZfK2nOmR5G0o3RGWtyDKzM33DxM1poA5MWLdp7oaoQLocWZLl408q19G7j+LEhRNQKLaH26V6G23inZCnPe9t9cCrPJNlfhNfeVU1+AOCkGsbnoVKBmVZK7KtlrgJknYKib43Ol3WLlkx6FFRsz3Tsixp5Zq9axanXzztfG/DOyFPe97b6oFXaTDK/CZavujgowjw0EN2nvHKGlvGivTVF540FnD8+I2tnrHpsnYp4md2aaZlUdLK5Utvwzshz3re2+iBV2kwyvwmWr74yqtNrd6YRlrQSSMa66svPN7KnT2bbmlUnS7rkLVeVKRdmWlZNHwwq1w+9DZEtZVlT65jampK19bWKv/eofu6EZrYnaxKHvKCTqxiupJMjUanWeRl0oj+xiNr3TZFp/ObPk+biMgFVb1hzo53Fjng/9hXHtEB3JmZ/GCPJogaiZubgaG5sOBIvRfdwq2oONcRy6SexIED5rpsPo8dtEib0/ld6W1E8VLIu45pI6yKsRm6eDY3A1/9s88C585l5KXNblLSQEI8rKfIegNJ6VURy7hlUcXSSPuNr2MHLWN6Or8vMzpDKOQOYtJlWrVRCI3EhYVAxIfDDI1ru/ufNZAwM1NuvYEyqzDawNexg5YxOZ3fpxmdIRRyB8kywspqZp1GYXo6EPJz53IMwra7/3FxC/MQFmJpqdh6A3EXTd4qjLao6Ev0zaqsg8np/K5GpmRBIXeQJCMsFNz19XKaGW8UJievbQ4N5DcKWQZhmKcHJmdwoKnuf1pLk+TOiBa6yHoD0QYIKLYKo0tktMI+WpV1MeW7LuOScaWxpJAbxKSbOKpTccOxjBcgKsSTk9cHVERXk81qFJIMwmiePjExjeeXV3BgY9Vs97+uXyipJYzmz2f/c07dNG1V1hEwV8QvjaIuGZcaS2+FPE00yx43lY+4SJp0E8cNx4cfDozHMkERoceh6mqyWXna2gJ+d2MaB+YN38RVXTZJMy2TRM9n/3NO3TS5SFYoYJvbm+j1enjy/icxd3Cu0G9PXjiJY08fw47uYHd/t7M9hSLWvUsuGC+FPO25LHvcZD5EArds5qAg6kWQRFeBrZL/pHSKriabl1YjxmyVkyRd7CzR8zWWNT7gu74eROqMgv2np5ubXbl6cRWb25sYYojhcIhjTx/DgVsPFFrb+9GnH8X2MBiE3tzZ9ML/nEbSYmFL55as9DS8FPK057Ls8SIU3by917sWQBENnihiGOZhynBMS6dKerXyVLQ1q3KSpIvtqwslbyR6ZSVYVe3UqUDEh8PgRhyFW05PNxPvPLN/Br1eD8PRgPL2cLvQEq+rF1ev/gYA+tL3wv+cRtQFM3nTJI4/c9yem0VVW/87ePCg1uH8edU9e1T7/eD1/Plqx6ueJ+3zEydUFxeD40m/XVwM3gPB6+JitbKH5yjzmTNUvRh10/et4orWU/SmCv+q3lwlOLF2QvuP9xULUCxAd39it55YO6GLzy3q+fXkvJ5fP697ntijvYWe7vrnu/TE2onU9MPv9h/v654n9qSmmfX7rLyYZvG5xav10X+8r4vPNVP/ANY0QVO9tMizLMsyx/PIs+Sz0o37o00Yhp2Yyd10qGLWTZDm63Kx4orWU3zmVq9X6uaqavXOHZzDVy5/5eoSrz/c+WHg+x7upPrNy8R1x/3PZ148UzifTQxC5tVTk2MSRfBSyIH057Ls8ZAqQQ1ZPd+k39Z1kWQ9203rYxqNr1NeJQNl/N62Ki6PovUUD0sqsSBOXcGLLvEqItgebkOhN/jN4yJY5BxRYez3+jj1wilsD7cL5dP0IGSRerK94qO3Qm6KweCamzG+XEdeDHWWIVfWMCxC1rNdRh9NRfC0sk658QzEcNV/XqaeKt5UdQUv7iOODmTu6M7VTRhCEez3+jh699FCu+xE017/wTqe+vJThfNZ1zqONzxF68nmGizeCbnJMMJQB8IdeIDiQQ1FDDnTARFZz3bR596kJ6GyMWuqYkxY0y6HIDYcUWPCHRAXr2hoYXz97p2dHZy4cAKnXzydaNUmWe6hRR/d3CEvn3Ws4yTr27bbpAheCblpd2aoA6GIRyNO8rBlyGU920Wee5OeBOvGbFIGqrT0voYg1sS0O2Du4BwO3HrghvQm+hN4c/tN6OhfklWb5b6oks+q1nGS9T3/gXnnN8rwSsizRKhufHa/Dxw9WjxG22VDLoui4lukPq3XQTwDgJsDlw6RZvWaIp5eKMJnXjxznZ87btXmuS/acltk7RTkooBfJSmUpem/quGHpsMLw9/mRZ65GJ1Wh7zyGIkQtFFpJuI7O0zdkD4T508LCaySt6ZCDNsOXSwDmg4/FJE+gDUA31XVB0ylGyXNAoxb6mfOFLcSi0SzdM3IyytzbfeLrUqz7utxhJTulO0p5VlWbVn3SZPrnDhvfSdg0rXyUQAvAXi7wTRvIEmE4i6SpAiUkLIuGFej05qkth7aqjRTvp6mF+YxmW6J6cOub7RQRkCbapRcn02aSpKZXvYPwLsArAD4KwB+N+/7dWd2JhH25B95JL13XcVl0PRERFep5RnxudKaynsT6eZNHxYJHojoTyq4DWy7ZNrKU5E0bbtd0LBrZRnAPwbwtrQviMgcgDkA2Ldvn6HTXiO01AcD4PTpZGsyaii++WbgginiemlqQM/lTaRrBXJYHwWtQZ3eRNGFeUz1UpLSnJkJuqU7O0E41qlT143gp1m9SRsXh+9tu2SSaGICTl45XVq2Nk5tIReRBwC8rqoXRGQm7XuqehLASQCYmprSuudNI0tDcu7xzDRNa1EXfe/X4WtIX1m/UtF1jJvw36dNHz56FDhxIrjJt7dzG424QC0fXr5uAajlw8tOxlGb9mXnuZ6iQr+5vYmF1QUszCw4IeYmLPL3A/jrInI/gB8B8HYR+Q1V/UUDaVcia5p+yXu8MUwZaC5b9V5SZmZVdEpw3jrGTfRS0tKcnU3vliYQt0TPfvPsde833thoPY7ahq86z8oPhT5cwvfZ7zyLc+vnSlnmTZVLVM0ZxyOL/Jc1J2plampK19bWjJ23DK5Ywiby4UpZxo6kKcHhOsbDoRsXo0QLn2eRN+VCSBM1Uy6MJkRzcGmAhdUFPPudZzHUIfrSx8P3PIx9P7Yv9zwmyiUiF1R1Kn7cqwlBJnDFfWsiH+MYUQPAfjckaUrw7t3A8nKpRasapYRrK8kSTZqhaZIsUTPhk2/Knz29dxoLMws4t36u9IJeTY41GBVyVV0FsGoyzSZwxX1bNx/eh01XEWQXuiF1pgS3TFGrNGlGZvi+Ccs2S9RMrG3SpGhWXdCryTVbxs4i7xKu9C4qUVWQXeiGeFLxJqzSptwcWaKW5asu2qg0vdBVlQW9mlzqttNC3sb8iybPVQRXehelqSrIrnRDPKh4E1Zpk26OLFFLikgps3FzW+uDlz1PU7NGOyvkTfTAm9jc2ba71xpVBdkTa7gQSRff4A1hwipt0s1RRtSqbNzc1lR7F6b0d1bI25p/MT1d/VymGhsvG4M6guyBNZxL0sUHjFofJqxSE2mYagzKbNwcpc1QxvBckzdNYuONjdbCJzsr5GUMvqJCmJZmVePSRGPjwthfZbogyFVJuvhAsRuiRMttwlqsm4apxmD3rt3Y3N5Er9fDp+7/VKF02pyNGZ4rjDPvSS/XBWSKzgp5mXkdRYUwLc2qxmVdd+9gACwsXNt315kQRC+7CC2TdvGj0TDr60FdxgdjPGm5q+zVmUbVxqDN5QXCcw0R9ByGOmxtSYPOCjlQzOAraxVnzRot+zzV8S6Ez3PFzdObwyOhsUraxV9ZuTZj9Kmnghma0TpsOWqnqluiCUu4SmPQ5jZt8ZmfPem1tqRBp4W8CLaDIKp6F8LnORTx++4LrHPrmulCeKAvRMU5fB8OumxvJ9dhizdsHTF2ZaGttqJX4ueij7xlfA2CiD/PTog4YL9l9Im03ktWHbZ4w9YRY5c2LM6z5E0OhtqKYBl7IQf8HHNztgFyNmMFadO/n9Z7yavDlm7YOmLcpiVclCTBdnlp2jJQyFugKW1wtgFyNmM5tO3fz7O8LddhVTGOCub8B+YbzWNRazpNsF1xAdWFQt4wHPtzjLY3f8jCg95LWVeBjXC/IudKE2yXXEB1oJA3DMf+HCIa6tPrAU8+CczNXfvchn/fAcvbJDbC/eosWOWiC6gKFPKG4difQ6yuXovXHA6BY8eAAwfqTwggV7ER7ld3waqmByjbmFlqdGOJotjcWMIGnB/jCIMBcO+9QWgfEFjlTzwBzDfrxy1MR24UG1PiXbWmTbuauLGERTrWe/aX6enAnXLsWODr2r272S5SGWHu0GCKaQs3S6xdWLAqi7ZcTRRyMl7MzQXulKYt36ylMpPOzcGURHwPD2zL1dRZIe9IL/U6XCiTC3moTRtdpLRFsdKsbg6mJOJyeGCeWyf8fPnwcuOzPDsp5F1cH9yFnrcLeWidqjdEkjBnWd0caE3E1fDAvJ5C2z2JTgp51V6qy0LlQs/bhTy0Sp0bIk2Ys6xuDqbcgKvhgXk9hbZ7Ep0UcpvrgzeFCz3vxvLgajeo7g0RF2ZHrO40l4CrESBNDGjWLWteT6HtnkRnww993aA9i6plqqMb8d9XSi/rR/FKX14GNjbcEHXXb4gKpHX5fR9ULENTG0qX/bwKYxd+2Pb64G1Qtkx1dSjt96XqJS8TUat3czMIDRwO3RBO12+ICqR1+V0eVDSNqbLm9RTaDI3stXIWj5ievjY/ZGkp0CFfSQucaOv3hRIJ/TX9fjBBZ2en5gkNE94QHRBx4FqXvy/967r8acd9YHBpgKVzSxhcGiS+j5NU1rzfuE5nLfIkiroFutKjruvTNuITz0skavVOTgLHjzMEr0HSBg9dHVTMI+4mWT68jOPPHM90m8TLCsB7t1JtIReRvQDOAPgpAEMAJ1X1k3XTNU0ZcXZ50LMMdT0DRjwLRRKJ+mvamKwz5qR1+V2fJRkl9D+v/2D9OjfJ2W+eLeQ2iZZ16dyS924lExb5NoB/pKpfFpG3AbggIl9U1W8aSNsYZcTZhQgRU9SNaDMSEVcmEYbgkRyiVni/18eu3i5gCEz0J/Dzd/w8zq2fKxUt4mqsehlqC7mqXgZwefT/PxWRlwDcBsApIS8jzh0c4yKkM0QHKzEEHr7nYez7sX1XXUIHbj1QykXkq1spitHwQxHZD+A5AHeq6p/EPpsDMAcA+/btO/jqq68aO29RXA1XJoQEFAnZqxI+6GqMfFnSwg+NCbmIvBXAfwPwL1T1d7K+O27L2BJC8ikj0GWEuUsx8mlCbiT8UETeAuAsgM/kiTgJGAz8D28kfuBLaF1SfHca03unMf+B+UKCXCZdXzERtSIAfg3AS6r6q/Wz1H2yImjo/iEm8ckabWrQsQuDmXmYiFp5P4C/A+BrIvLC6Ng/VdWnDaTdSdIiaLoSv07cwacZm00NOnZhMDMPE1Er/x2AGMjL2JAWQdOV+HXiDr5Zo03FsvsUI1+FsZrZ6Qpp4Y1dil8nbjAO1ijp8OqHvkIfOSFu4kII49itfugrnNjoEWx1xwbXB40p5IRUgSPTY4Xrg8ZcxpaQKhhZ49dtfIk/bwPXl/mlRT6m0CtQExsj0y1eNNddCW3j+qAxhXwMoVfAAG2vrNbyRXPdlWADl0MYKeRjCOPVDdHmyHTLF823+PNxh0I+hjBe3UNavmiuuxLI9TCOfEyhj9xDeNHGnsaXsS2DSSHnvU0IGRc6OSGIg3aEEOJ5HPkYhPISQkguXgt5OP7T73PQjjiGBzuHcMJPd/DatcJNkomTeODz44SfbuG1RQ4Ez8f8vHPPCRlXBgNgYQHY3HTa5zcO25+NE15b5KQejPgxTGiJb24CwyHQ6znr8+OEn25BIR9TPOj9+0c4+h6K+H33Bda5gxXLCT/dgkI+pnCafgPEZ186KuIhLq8dQspBIR9TOE2/ATj6TixBIR9TqDkNwS2eiAUo5GMMNYeQbuB9+CEhhIw7FHJCCPEcI0IuIodF5GUReUVEHjORJvljXxoAAAaLSURBVCGEkGLUFnIR6QN4EsCHANwB4BdE5I666RJCCCmGCYv8vQBeUdVvq+oWgM8CeNBAuoQQQgpgQshvA3Ap8v610bHrEJE5EVkTkbUrV64YOC0hhBDAjJBLwrEbth1S1ZOqOqWqU7fccouB0xJCCAHMCPlrAPZG3r8LwPcMpEsIIaQAJoT8DwC8R0RuF5EJAB8G8HkD6RJCCClA7ZmdqrotIscA/FcAfQCfVtVv1M4ZIYSQQhiZoq+qTwN42kRaVeC62oSQccb7tVa4rjYhZNzxfop+0rrahBAyTngv5OG62v0+19UmhIwn3rtWuK42IWTc8V7IAa6rTQgZb7x3rRBCyLhDISeEEM+hkBNCiOdQyAkhxHMo5IQQ4jkUckII8RwKOSGEeA6FnBBCPIdCTgghnkMhJ4QQz6GQE0KI51DICSHEcyjkhBDiORRyQgjxHAo5IYR4DoWcEHIdg0sDLJ1bwuDSwHZWSEE6sbEEIcQMg0sDHDpzCFs7W5joT2BldgXTe7lri+vQIieEXGX14iq2drawozvY2tnC6sVV21kiBaCQE0KuMrN/BhP9CfSlj4n+BGb2z9jOEikAXSuEkKtM753GyuwKVi+uYmb/DN0qnlBLyEXk3wD4awC2APwhgIdU9f+YyBghxA7Te6cp4J5R17XyRQB3qurPAPifAObrZ4kQQkgZagm5qn5BVbdHb38PwLvqZ4kQQkgZTA52HgXwX9I+FJE5EVkTkbUrV64YPC0hhIw3uT5yEXkWwE8lfPQxVf3c6DsfA7AN4DNp6ajqSQAnAWBqakor5ZYQQsgN5Aq5qt6X9bmIHAHwAIBDqkqBJoSQlqkbtXIYwD8B8JdV9Q0zWSKEEFIGqWNEi8grAHYD2Bgd+j1VfaTA764AeLXiaW8G8EcVf+saLIubsCxuwrIAf15Vb4kfrCXkNhCRNVWdsp0PE7AsbsKyuAnLkg6n6BNCiOdQyAkhxHN8FPKTtjNgEJbFTVgWN2FZUvDOR04IIeR6fLTICSGERKCQE0KI53gl5CJyWEReFpFXROQx2/kpi4hcFJGvicgLIrI2OvYOEfmiiHxr9PoTtvOZhIh8WkReF5GvR46l5l1E5kfX6WUR+at2cn0jKeVYEJHvjq7LCyJyf+QzJ8sBACKyV0S+JCIvicg3ROSjo+M+Xpe0snh3bUTkR0Tk90XkxVFZHh8db+66qKoXfwD6CNY8fzeACQAvArjDdr5KluEigJtjx/41gMdG/38MwL+ync+UvN8L4B4AX8/LO4A7RtdnN4DbR9etb7sMGeVYAPDLCd91thyj/L0TwD2j/78NwVLSd3h6XdLK4t21ASAA3jr6/1sAPA/gfU1eF58s8vcCeEVVv62qWwA+C+BBy3kywYMATo/+fxrAz1nMSyqq+hyAP44dTsv7gwA+q6qbqvodAK8guH7WSSlHGs6WAwBU9bKqfnn0/z8F8BKA2+DndUkrSxoul0VV9f+O3r5l9Kdo8Lr4JOS3AbgUef8asi+0iyiAL4jIBRGZGx37SVW9DAQ3M4BbreWuPGl59/FaHRORr45cL2GX15tyiMh+AH8RgfXn9XWJlQXw8NqISF9EXgDwOoAvqmqj18UnIZeEY77FTr5fVe8B8CEAj4rIvbYz1BC+Xat/D+AvALgbwGUA/2503ItyiMhbAZwFcFxV/yTrqwnHnCpPQlm8vDaquqOqdyPYbOe9InJnxtdrl8UnIX8NwN7I+3cB+J6lvFRCVb83en0dwH9G0H36voi8EwBGr6/by2Fp0vLu1bVS1e+PHrwhgKdwrVvrfDlE5C0IhO8zqvo7o8NeXpeksvh8bQBAgz2MVwEcRoPXxSch/wMA7xGR20VkAsCHAXzecp4KIyI/KiJvC/8P4IMAvo6gDEdGXzsC4HN2cliJtLx/HsCHRWS3iNwO4D0Aft9C/goRPlwj/gaC6wI4Xg4REQC/BuAlVf3VyEfeXZe0svh4bUTkFhH58dH/9wC4D8D/QJPXxfYIb8nR4PsRjGb/IYIdiqznqUTe341gZPpFAN8I8w9gEsAKgG+NXt9hO68p+f9NBF3bHyKwID6SlXcAHxtdp5cBfMh2/nPK8R8BfA3AV0cP1TtdL8cob38JQRf8qwBeGP3d7+l1SSuLd9cGwM8A+Mooz18H8M9Gxxu7LpyiTwghnuOTa4UQQkgCFHJCCPEcCjkhhHgOhZwQQjyHQk4IIZ5DISeEEM+hkBNCiOf8f4haQX6BdMXGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "std2 = 2\n",
    "data2 = []\n",
    "for ii in range(3):\n",
    "    data2.append(stats.norm(centers[ii], std2).rvs(100))\n",
    "    plt.plot(np.arange(len(data1[ii]))+ii*len(data2[0]), data2[ii], '.', color=colors[ii])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Note:** In both cases the means have the same difference, but the variance is much larger in data2!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ANOVA with Sample Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Get and sort sample data\n",
    "\n",
    "*Twenty-two patients undergoing cardiac bypass surgery were randomized to one of three ventilation groups:*\n",
    "    \n",
    "  - *Group I: Patients received 50% nitrous oxide and 50% oxygen mixture continuously for 24 h.*\n",
    "  - *Group II: Patients received a 50% nitrous oxide and 50% oxygen mixture only dirng the operation.*\n",
    "  - *Group III: Patients received no nitrous oxide but received 35-50% oxygen for 24 h.*\n",
    "    \n",
    "*The data show red cell folate levels for the three groups after 24h' ventilation.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Get the data\n",
    "inFile = 'altman_910.txt'\n",
    "url_base = 'https://raw.githubusercontent.com/thomas-haslwanter/statsintro_python/master/ipynb/Data/data_altman/'\n",
    "\n",
    "url = url_base + inFile\n",
    "data = np.genfromtxt(urllib.request.urlopen(url), delimiter=',')\n",
    "\n",
    "# Sort them into groups, according to column 1\n",
    "group1 = data[data[:,1]==1,0]\n",
    "group2 = data[data[:,1]==2,0]\n",
    "group3 = data[data[:,1]==3,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Levene-test\n",
    "A Levene-test and/or a normality test should be made before applying a oneway ANOVA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: the p-value of the Levene test is <0.05: p=0.045846812634186246\n"
     ]
    }
   ],
   "source": [
    "# check if the variances are equal with the \"Levene\"-test\n",
    "(W,p) = stats.levene(group1, group2, group3)\n",
    "if p<0.05:\n",
    "    print('Warning: the p-value of the Levene test is <0.05: p={0}'.format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### One-way ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The results from the one-way ANOVA, with the data from Altman 910: F=3.7, p=0.04359\n",
      "One of the groups is significantly different.\n"
     ]
    }
   ],
   "source": [
    "F_statistic, pVal = stats.f_oneway(group1, group2, group3)\n",
    "\n",
    "print('The results from the one-way ANOVA, with the data from Altman 910: F={0:.1f}, p={1:.5f}'.format(F_statistic, pVal))\n",
    "if pVal < 0.05:\n",
    "    print('One of the groups is significantly different.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Elegant alternative implementation, with pandas & statsmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                df        sum_sq      mean_sq         F    PR(>F)\n",
      "C(treatment)   2.0  15515.766414  7757.883207  3.711336  0.043589\n",
      "Residual      19.0  39716.097222  2090.320906       NaN       NaN\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame(data, columns=['value', 'treatment'])    \n",
    "\n",
    "# the \"C\" indicates categorical data\n",
    "model = ols('value ~ C(treatment)', df).fit()\n",
    "\n",
    "print(anova_lm(model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
